Force chains and contact network topology in packings of elongated particles 
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By means of contact dynamic simulations, we investigate the contact network topology and force 
chains in two-dimensional packings of elongated particles modeled by rounded-cap rectangles. The 
morphology of large packings of elongated particles in quasistatic equilibrium is complex due to 
the combined effects of local nematic ordering of the particles and orientations of contacts between 
particles. We show that particle elongation affects force distributions and force/fabric anisotropy 
via various local structures allowed by steric exclusions and the requirement of force balance. As 
a result, the force distributions become increasingly broader as particles become more elongated. 
Interestingly, the weak force network transforms from a passive stabilizing agent with respect to 
strong force chains to an active force-transmitting network for the whole system. The strongest 
force chains are carried by side/side contacts oriented along the principal stress direction. 

PACS numbers: 45.70.-n,83.80.Fg,61.43.-j 



I. INTRODUCTION 

Most remarkable properties of granular materials are 
closely related to their specific disorder induced essen- 
tially by steric exclusions and the force balance condition 
for each particle. The broad and strongly inhomogeneous 
distribution of contact forces, as a hallmark of granular 
disorder, has been a subject of extensive investigation ^ 
[To| . In spite of particle mobility and disorder, granular 
materials exhibit a finite shear strength due to a gen- 
uine anisotropic two-phase organization of the contact 
network involving strong force chains propped by weak 
forces pTli- 

The robustness of these micro-structural features with 
respect to particle geometry and interactions has been 
addressed only recently by discrete-element numerical 
simulations. For example, it is found that in highly poly- 
disperse systems the force chains are mainly captured 
by large particles so that the shear strength of a non- 
cohesive granular material is practically independent of 
particle size distribution 12|. As another important ex- 
ample, a parametric study shows that when the particles 
interact by both sliding friction and high rolling resis- 
tance at their contacts, the nature of the weak network 
is affected by the formation of columnar structures which 
do not need to be propped by a particular class of weak 
contacts p^ . 

The particle shape is another major characteristic of 
granular material. Most applications of real granular ma- 
terials involve some degree of deviation with respect to 
simple circular or spherical shapes often used in simula- 
tions by the discrete-element method. While the numeri- 
cal treatment of large packings of complex particle shapes 
was until very recently out of reach due to demanding 
computational resources, there is presently considerable 
scope for the numerical investigation of complex granular 



packings. This is not only due to the enhanced computer 
power and memory but also because during more than 
two decades of research in this field, many properties of 
granular media have been investigated for model packings 
composed of circular and spherical particle shapes. Such 
properties provide thus a rich guideline for the analysis 
of specific behaviors arising from particle geometry. 

Systematic studies of particle shape dependence in 
granular materials have been recently reported for polyg- 



onal/polyhedral [14 - 
non-convex shapes 
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elliptical/ellipsoidal 19-22] and 
23[ . The force chains are found to be 
reinforced in packings of polygonal and polyhedral parti- 
cles leading to enhanced shear strength [l^ [13 [131 • The 
effect of shape elongation was investigated for packings 
of rectangular-shaped particles deposited under gravity 
[25| . The preparation under gravity has strong influ- 
ence on the particle orientations and thus on the force 
distributions. On the other hand, a systematic study of 
the shear behavior of 2D packings of rounded-cap rectan- 
gles (RCR) under homogeneous boundary conditions in- 
dicates that the shear strength increases with elongation 
whereas the packing fraction varies unmonotonically 23) , 
as also found for packings of ellipsoidal shapes [20, [2l| . 
In all reported cases, the networks resulting from various 
shapes appear to be highly complex and hardly amenable 
to simple statistical modeling. 

In this paper, we use contact dynamics simulations 
to investigate the contact and force networks in sheared 
granular packings of RCR particles with increasing as- 
pect ratio in 2D. We focus more specifically on the orga- 
nization of the contact force network in correlation with 
the fabric anisotropy described in terms of branch vec- 
tors joining particle centers. Our data reveal a bimodal 
force network as in disk packings but with qualitatively 
different roles of fabric and force anisotropics. This be- 
havior involves a short-range nematic ordering of the par- 
ticles with side/side contacts that capture stronger force 
chains. On the other hand, the friction mobilization is 
shown to be anisotropic and it plays a major role for the 
■fltatttiatgf I of elongated particles. 
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In the following, we first briefly describe the numeri- 
cal procedures, which are essentially the same as those 
reported in [22|. Then, we analyze the branch vectors 
and their correlations with the contact forces. Finally, 
we present a detailed analysis of the partial stresses and 
fabric anisotropics sustained by force sub-networks. We 
conclude with salient results of this work its possible 
prospectives. 
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FIG. 1: Shape of a Rounded-Cap Rectangle (RCR). 



II. MODEL DESCRIPTION AND NUMERICAL 
SIMULATIONS 

The simulations were carried out by means of the con- 
tact dynamics (CD) method with irregular polyhedral 
particles. The CD method is a discrete element approach 
for the simulation of nonsmooth granular dynamics with 
contact laws expressing mutual exclusion and dry fric- 
tion between particles without elastic or viscous regular- 
ization [10, 26 32]. Hence, this method is particularly 
adapted for the simulation of perfectly rigid particles. 
Nonsmoothness refers to various degrees of discontinu- 
ity in velocities arising in a system of rigid particles. In 
this method, the equations of motion for each particle 
arc formulated as differential inclusions in which velocity 
jumps replace accelerations [1^. The unilateral contact 
interactions and Coulomb friction law are treated as com- 
plementarity relations or set-valued contact laws. The 
time-stepping scheme is implicit but requires explicit de- 
termination of the contact network. Due to implicit time 
integration, inherent in the CD method, this scheme is 
unconditionally stable. 

At a given step of evolution, all kinematic constraints 
implied by lasting contacts and the possible rolling of 
some particles over others are simultaneously taken into 
account, together with the equations of dynamics, in or- 
der to determine all velocities and contact forces in the 
system. This problem is solved by an iterative process 
pertaining to the non-linear Gauss-Seidel method which 
consists of solving a single contact problem, with other 
contact forces being treated as known, and iteratively up- 
dating the forces and velocities until a convergence crite- 
rion is fulfilled. The iterations in a time step are stopped 
when the calculated contact forces are stable with respect 
to the update procedure. To check convergence we thus 
use the relative variation of the mean contact force be- 
tween two successive iterations. We require this relative 
variation to be below a given value which sets the preci- 
sion of the calculation. In this process, no distinction is 
made between smooth evolution of a system of rigid par- 
ticles during one time step and nonsmooth evolutions in 
time due to collisions or dry friction effects. The unique- 
ness of the solution at each time step is not guaranteed 
by CD method for perfectly rigid particles. However, by 
initializing each step of calculation with the forces calcu- 
lated in the preceding step, the set of accessible solutions 
shrinks to fluctuations which are basically below the nu- 
merical resolution. In this way, the solution remains close 



to the present state of forces. 

For our simulations, we used the LMGC90 which is 
a multipurpose software developed in Montpellier, ca- 
pable of modeling a collection of deformable or unde- 
formable particles of various shapes (spherical, polyhe- 
dral, or polygonal) by different algorithms [13, [13 • 

A. Simulation of RCR particles 

We model the RCR particle as a juxtaposition of two 
half-disks of radius R' with one rectangle of length L 
and width 2R'\ see Fig. [TJ The shape of a RCR parti- 
cle is a circle of radius R' for L = 0. The aspect ratio 
a — [L -\- 2i?')/(2i?') is 1 in this limit and increases with 
L for a fixed value of R' . In this paper, we use an alter- 
native parameter describing the deviation of the particle 
shape from a circle. Let R be the radius of the circle 
circumscribing the particle. We have R = L/2 + R' . The 
radius R' is also that of the inscribed circle. Hence, the 
deviation from a circular shape can be characterized by 
Ai? — R — R' — L/2. We use the dimensionless parame- 
ter 7] defined by 



It varies from 77 = 0, for a circle, to 1 corresponding to a 
line. We will refer to 77 as the elongation parameter as in 
rock mechanics (ssj . 

The contacts between RCR particles belong to differ- 
ent categories, namely cap-to-cap (cc), cap-to-side (cs) 
and side-to-side (ss); see Fig. [5] Side-to-side contacts re- 
sults from contacts between two rectangles as well as two 
contacts resulting from cap-to-side. In the CD method 
the case of side-to-side contacts for rectangular particle is 
represented by two points. Hence, for RCR particles, ss 
contact is composed of four point contacts : two points 
due to the rectangle-rectange interface and two points 
due to the cs contacts. In the iterative procedure of 
determination of the contact forces and velocities, the 
points representing the contact between two particles are 
treated as independent points but the resultant of the 
calculated forces are attributed to the contact with its 
application point located on the contact plane. 

The detection of line contacts between rectangles 
was implemented through the so-called shadow overlap 
method devised initially by Moreau [12, 113| for polygons. 



FIG. 2: Representation of cap-to-cap, cap-to-side and side- 
to-side contact and they will be referred as cc contacts, cs 
contacts and ss contact, respectiveley. 



The reliability and robustness of this method have been 
tested in several years of previous applications to granu- 
lar materials [15, 16, 32. 34-37]. This detection procedure 
is fairly rapid and allows us to simulate large samples 
composed of RCR particles. 



B. Packing preparation and bi-axial test 

We prepared 8 different packings of 13000 RCR par- 
ticles with rj varying from to 0.7 by steps of 0.1. The 
radius R of the circumscribing circle defines the size of 
a RCR particle. In order to avoid long-range ordering in 
the limit of small values of rj, we introduce a size poly- 
dispersity by taking R in the range [Rrnin,Rmax] with 
Rmax — '2Rmin with a Uniform distribution in particle 
volume fractions. 

All samples are prepared according to the same pro- 
tocol. A dense packing composed of disks (77 = 0) is 
first constructed by means of a layer- by- layer deposition 
model based on simple geometrical rules p38l - t40| . The 
particles are deposited sequentially on a substrate. Each 
new particle is placed at the lowest possible position at 
the free surface as a function of its diameter. This proce- 
dure leads to a random close packing in which each parti- 
cle is supported by two underlying particles and supports 
one or two other particles. For 77 > 0, the same packing 
is used with each disk serving as the circumscribing circle 
of a RCR particle. The RCR particle is inscribed with 
the given value of 77 and random orientation in the disk. 

Following this geometrical process, the packing is com- 
pacted by isotropic compression inside a rectangular 
frame of dimensions Iq x /iq in which the left and bottom 
walls are fixed, and the right and top walls are subjected 
to a compressive stress ctq. The gravity g and friction 
coefficients /i between particles and with the walls are 
set to zero during the compression in order to avoid force 
gradients and obtain isotropic dense packings. Fig. [3] 
displays snapshots of the packings for several values of rj 
at the end of isotropic compaction. 

The isotropic samples were sheared by applying a 
downward displacement on the top wall at constant ve- 
locity for a constant confining stress acting on the lat- 
eral walls; see http://cgp-gateway.org/ref010 for video 
samples. During shear, the friction coefficient fi between 




FIG. 3: Examples of the generated packings at the initial 
state. 

particles was set to 0.5 and to zero with the walls. The 
strain rate was low so that the shearing is basically of 
quasi-static nature. The internal angle of friction ip at 
every stage of shearing is given by 



with CTi > (72 are the principals values of the stress tensor 
(T. sin(y9 increases with shear strain and saturates to a 
constant value corresponding to the critical state which 
is a state independent of the initial configuration of the 
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FIG. 4: Local contact frame (n, t) 



packing. The critical-state value of siniyj represents the 
shear strength of the packing and it increases linearly 
with 7j from - 0.3 (for 77 = 0.0) to - 0.51 (for 77 = 0.7) 

In granular media, the expression of stress tensor cr in 
the volume V is an arithmetic mean involving the branch 
vectors (joining the centers of the two touching parti- 
cles) and contact force vectors at contact c, and it is 
given bv (illi^: 



(3) 



For the analysis of stress transmission from a particle- 
scale viewpoint we need a statistical description of these 
quantities. 

A common approach used by various authors is to ex- 
press branch vectors and contact force orientations in 
terms of the contact direction, i.e. in the local contact 
frame (n,t), where n is the unit vector perpendicular 
to the contact plane, and t is an orthonormal unit vec- 
tor oriented along the tangential force; see figure SI The 
components of the branch vector and contact force are 
expressed in the following frame: 



£ = inTl + itt, 

f = fnn + ftt, 



(4) 



where and it are the normal and tangential compo- 
nents of the branch vectors, and /„ and ft the normal 
and tangential components of the contact force. Remark 
that only for disks or spherical particles we have £ = in 
where i is the length of the branch vector. 

In the following we study the shapes of the normal and 
tangential force and branch distributions in the residual 
state. 



III. DISTRIBUTIONS OF CONTACT FORCES 
AND BRANCH VECTORS 

A specific feature of the contact network of a packing of 
elongated particles is that the length £ of branch vectors 



strongly varies throughout the network depending on the 
relative particle orientations. From the definition of rj 
(Eq. dl])) and for given values of Rmin and Rmax, it is 
easy to see that 



Rn 



2#^(l-^),2 



(5) 



In our simulations, since Rmin /Rmax — 0.5, we have 
{l — ri)Rmax < i < 2Rmax- With increasing elongation 77, 
the range of i becomes significant and its statistics can 
be used as a meaningful characterization of the texture 
as a function of rj. On the other hand, the correlation of i 
with the total reaction force / between neighboring par- 
ticles seems to be a good descriptor of the organization 
of forces for particles of non circular shape. The branch 
vectors are also important as they enter the expression 
of the stress tensor given by Eq. ([3]). In Ref. [22|, a 
different point of view was adopted: the contact forces 
were projected along and perpendicular to the branch 
vectors and their statistics were investigated. The same 
framework was used for the decomposition of the total 
stress tensor. Here, we focus on the distribution of con- 
tact forces and their correlation with the branch vector 
as 77 is increased. 



A. Contact forces and friction mobilization 

The probability density function (pdf ) of normal forces 
normalized by the mean normal force (/„) is shown in 
Fig. [5] in log-linear and log-log scales at large strains 
(the data are cumulated from several snapshots in the 
critical state) for all simulated values of 77. As usually 
observed, in all packings the number of forces above the 
mean (/„) falls off exponentially whereas the number of 
forces below the mean varies as a power-law: 

e-.(.)(/./(/.» , /„>(/„), 

Pifn)^i ( f \Mv) ^ ,^ , (6) 




fn < {fn)j 



where ari(77) and Pniji) whose variations are shown in 
the insets as a function of rj. We see that a„ decreases 
with increasing 77, implying that the inhomogeneity of 
normal forces becomes higher as the particles become 
more elongated. On the other hand, /3„ declines from 
0.1 to —0.4 with 77 which means that the proportion of 
weak contacts (carrying a normal force below the mean) 
increases with elongation. The proportion of weak forces 
grows from 60% for 7; = to 70% for rj — 0.7. In other 
words, while the proportion of strong contacts declines 
with increasing 77, stronger force chains occur at the same 
time. 

Figure [6] shows the pdf P{ft) of tangential forces nor- 
malized by the mean tangential force (|/t|) in each pack- 
ing. These distributions show also an exponential falloff 
for the forces above the average force (|/t|) and a power 
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FIG. 5: Probability distribution function of normal forces /„ 
normalized by the average normal force (/„) in log- linear (a) 
and log-log (b) scales for different values of rj. 
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FIG. 6: Probability distribution function of tangential forces 
ft normalized by the average tangential force (ft) in log- linear 
(a) and log-log (b) scales for different values of i]. 



law for the forces below 

the corresponding exponents a* (77) and /3t(?7) decreasing 
with r/. We observe that, in contrast to q;„ and /3„, the ex- 
ponents at and /3t saturate beyond rj = 0.4. This means 
that the friction forces do not follow the normal forces as 
77 increases. In other words, the most mobilized (largest) 
friction forces do not occur necessarily at the contacts 
where the normal forces are higher. 

In order to investigate the properties of friction mo- 
bilization, we consider the friction mobilization index 
Im = |/t|/M/«- Its average Im = (^) increases from 
0.4 for 77 = to 0.6 for rj = 0.7 as we see in Fig. [T] This 
increase underlies to a large extent the increase of the 
shear strength with 77, as we shall see below in Sec. IIVI 
However, the friction force is not uniformly mobilized at 
all contacts. Fig. [5] shows a map of weak (/„ < (/„)) 
and strong (/„ > (/n)) normal forces, represented by 
the thickness of vectors joining the particle centers to 
the contact points, and the corresponding values of /,„, 
represented by circles of diameter proportional to /,„ for 
77 = 0.1 and77 = 0.7. Visual inspection reveals that most 
mobilized contacts belong to the weak force network. In 
fact, the average friction mobilization Imf defined as the 
average by force class, plotted as a function of /„ in Fig. 
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FIG. 7: Friction mobilization Im averaged in the steady state 
as function of r;. 



|9]for all values of 77, declines as /„ increases. We also see 
that the friction mobilization increases with 77 at all force 
levels. 

Figure [10] displays the pdf of /,„ for different values of 
77 in the critical state. For the disks, the pdf is a nearly 
decreasing linear function of /„, which means that the 
proportion of weakly mobilized contacts is larger than 
that of strongly mobilized contacts. As 77 is increased, the 
distribution becomes more uniform, and at even larger 77 
a class of highly mobilized contacts (with I„i close to 1) 
appears whereas the distribution is nearly uniform for all 
other contacts. This class belongs to weak force network 
as was shown previously, so that not only the friction 
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FIG. 8: A snapshot of the force-bearing particles at = 0.2(a) 
and 77 = 0.7(b) and normal forces represented by the thickness 
of the segments joining the particle centers to the application 
point of the force. The strong and weak forces are in back and 
red, respectively. The diameter of yellow circle is proportional 
to Im at the contact. 
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FIG. 9: Friction mobilization Imf 
class, as a function of fn for all rj. 



as the average by force 



mobilization /,„ but also the number of highly mobilized 
contacts are larger in the weak force network. A class of 
very weak forces was also evidenced in [43| in a packing 
of disks deposited under gravity and tilted towards its 
angle of stability. This subclass of the weak network 
can be defined as the class of contacts where the normal 
force is below the mean but the friction is highly or fully 
mobilized. 

This enhanced friction mobilization implies that the 
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FIG. 10: Probability distribution function of the friction mo- 
bilization index /,„. 




FIG. 11: A snapshot of the force-bearing particles at r; = 0.7 
and normal forces represented by the thickness of the seg- 
ments joining the particle centers to the application point of 
the force. The color level for the particles is proportional 
to the orientation of the major particle axis for the particles 
with at least one side/side contact. The particles having no 
side/side contacts are in gray. The strong and weak forces are 
in back and red, respectively. 



equilibration of the particles is more complex than in 
disk packings. In particular, the nematic ordering due 
to the "geometrical" chains of side / side contacts between 
particles means that the statistics of forces and the mobi- 
lization of friction are closely related to the equilibrium 
of such chains rather than that of individual particles. 
These chains are evidenced in Fig. [Tl]for 77 = 0.7 where 
the force bearing particles belonging to the chains are 
represented by a color level proportionally to their orien- 
tations. The friction needs to be highly mobilized inside 
the chains in order to ensure their stability. 



B. Branch vectors 

The branch vectors in a packing of elongated particles 
reflect both the relative orientations of the particles in 
contact and their size distribution. The latter may be in- 
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FIG. 13: Probability distribution function of the reduced 
branch- vector lengths Ir for all values of t) in the critical state. 



FIG. 14: A snapshot of cap-cap modes contact (blue), 
side-side modes contacts (red) and cap-side modes contacts 
(green) . 



tegrated out by simply dividing the branch vector length 
^ between two touching particles by the sum Ri + R2 of 
the radii i?i and R2 of their circumscribing circles. This 
reduced branch-vector length ^r = £/{Ri + R2) varies 
in the range [1 — ?7, 1]. We have = 1 at 77 = (for 
disks). For elongated particles, — ^ corresponds to 
a cap/cap contact between two aligned particles, Fig. 
n^ a). whereas Ir = 1 — corresponds to a centered 
side/side contact between two parallel particles. Fig. 
I12f b). Such contact configurations, when they exist, can 
be evidenced from the probability density function of £r 
and its possible modes at fj. = 1 or = 1 — V- 

Figure [T3] displays the pdf of reduced branch- vector 
lengths for all values of 77 in the critical state. These 
pdf's are nothing but normalized radial functions with 
£r varying in a limited range as only the touching parti- 
cles are considered. They are nearly similar for all values 
of r]. The first mode, centered on = ^ — V reveals 
the presence of a broad population of side/side contacts 
with a peak increasing in amplitude with 77 as displayed 
in Fig. [14] (b). We also observe a less pronounced mode, 
centered on — 1, corresponding to a distinct popu- 
lation of aligned cap/cap contacts, also marked in Fig. 
[13] (a). The intermediate mode occurs approximately at 
£r ~ 1— r//2 which is the midpoint of the interval [1— 77, 1]. 
This length corresponds to an orthogonal side/cap con- 
tact as shown in Fig. [T^ c). The presence of such a 
distinct mode, through decreasing in amplitude as 77 in- 
creases, is a clear proof of the occurrence of orthogonal 
layers some of which are observed in Fig. 1141 This mode 
is also characterized by a broad extension refiecting the 
intermediate angles between the orientations of touching 



particles. 

We expect the branch vectors lengths to be correlated 
with contact forces because of either the contact con- 
figurations they represent or simply the fact that force 
chains tend to be captured by larger particles (hence, 
longer branch vectors) [l^ • This correlation can be esti- 
mated with the Pearson coefficient, which for two random 
variables x and y is defined by the scalar 

^ {ix-{x))iy-{y))) 

^{{x-{xW)^{{y-{yWy 

Note that C = 1 corresponds to a full inter-dependence 
whereas C — Q means full statistical independence of the 
two variables. Fig. [15] shows the Pearson coefficients 
C/f^, between the force amplitude / and ir, as well as 
C/f, between / and I, as a function of 77. Both coeffi- 
cients decrease with 77 from positive values for 77 < 0.3 to 
negative values down to —0.22. The positive correlation 
(larger forces at longer branch vectors) is a consequence 
of the fact that the distribution of branch lengths at low 
values of 77 is governed by particles sizes. On the other 
hand, the negative correlation (larger forces at shorter 
branch vectors) refiects the effect of the increasing num- 
ber of side/side contacts as the particles become more 
elongated. 

Further insight into this force/branch- length correla- 
tion can be obtained from the average force amplitude 
calculated by taking the average force in a class 
of contacts in the interval [Ir — l^lr/'^^ir + as a 

function of tri shown in Fig. llGf a) for all values of rj. 
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FIG. 15: Correlation Cf^/ and Cf/ as a function of averaged 
in the critical state. 



This plot shows that for ah contact classes the associ- 
ated niean force {f)i^ is nearly equal to the global mean 
force (/) except to the class of the shortest branch vectors 
(side/side mode), which concentrates a mean force above 
(/), and the class of the longest branch vectors (cap/cap 
mode), which seem to carry a considerably lower force on 
the average. In his way, the rather weak correlation be- 
tween the reduced branch length and force appears here 
to be governed by the two afore-mentioned modes. In or- 
der to evidence the effect of particle size distribution, let 
us consider the average force amplitude {f)i as a function 
of i as shown in Fig. [TBl b) for all values of 77. For r] < 0.3, 
the contact force is on the average an increasing function 
of £. For disks {rj = 0), the variation of is a consequence 
only of the particle size distribution and, therefore, the 
increase of the mean force with £ means that the larger 
particles, involved in the longer branch vectors, capture 
higher forces. The same effect seems thus to under ly also 
the increasing mean force with rj for elongated particles 
with r] < 0.3. But at larger elongations, the trend is 
reversed and we see that the mean force declines as £ in- 
creases, reflecting thus the effect of the side/side contact 
mode as discussed previously. 



IV. WEAK AND STRONG FORCE NETWORKS 

The complex network of contact forces in a packing 
of elongated particles can also be analyzed by consider- 
ing the contribution of various classes of forces and/or 
branch vectors to stress transmission. Indeed, according 
to equation ^ , the stress tensor is expressed as an aver- 
age involving branch vectors and contact forces, so that 
partial summations allow one to define partial stress ten- 
sors that have been applied in the past to investigate the 
scale- up of local quanties For example, the subset of 
contacts carrying a force below a threshold, reveals the 
respective roles of weak and strong force chains with re- 
spect to the overall shear strength of granular materials 
In this section, we apply this methodology to ana- 
lyze the stress and other texture-dependent quantities in 
view of elucidating the effect of particle elongation. 




FIG. 16: Linear correlation between contact force / and 
branch length £ as a, function of rj. 



In what follows, we consider various fabric and force 
parameters for the "^-networks" defined as the subsets 
S{^) of contacts which carry a force below a cutoff force ^ 
normalized by the mean force (ie f/{fn) € [0,C])i where 
^ is varied from to the maximal force in the system. 
The weak network corresponds to 5(1) whereas the strong 
network is its complement. In section Hill we focused on 
scalar descriptors of granular texture such as the distri- 
butions and correlations of force magnitudes and branch 
lengths. Beyond these low-order quantities, the granular 
texture is characterized by a disordered but anisotropic 
structure of both the contact and force networks, which 
require higher-order description in terms of various fabric 
and force tensors. We analyze below different parameters 
pertaining to this tensorial organization of our packings 
as a function of ^ and for increasing elongation rj. 



A. Granular texture 

A relevant description of granular texture is given by 
the probability distribution P{n) of the contact normals 
n ; see Fig. |4l In two dimensions, the unit vector n 
is described by a single angle 6 e [0,7r]. The distribu- 
tion Pg {9) of contact orientations can be evaluated from 
the numerical data at different stages of its evolution. In 
our simulations, all numerical samples are prepared in 
an isotropic state so that Pg = I/tt in the initial state. 
This distribution evolves with shear strain and becomes 
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FIG. 17: Distributions of contact orientations (symbols) in 
polar coordinates for rj — 0.5 and several values of the force 
cutoff ^ together with their Fourier fits (|11[) (full lines). 

increasingly more anisotropic as the critical state is ap- 
proached. By restricting the data to those belonging to 
the ^-networks, we obtain a continuous family of distri- 
butions Pe{0,£) that describe the geometrical state of 
the system. In practice, however, such functions can be 
estimated with meaningful statistics only in the critical 
state where the data can be cumulated from independent 
configurations representing all the same state. 

Figure [T71 shows the distributions Pg{9,£) in polar co- 
ordinates for r] = 0.5 and for several values of ^. The 
distributions are similar with nearly the same privileged 
direction aligned with the principal stress direction 6^^ 
but with increasing anisotropy as a function of ^. They 
all can be approximated by their truncated Fourier ex- 
pansion [TTI.li3l.fi5j: 

Pe[e,0 = + COS 2(0-0,)}, (9) 

where adC) is the amplitude of contact anisotropy in the 
^-network. In practice, it is more convenient to estimate 
a^^) through the partial fabric tensors F{^) defined by 

Fo.m^- r n^{0,O^I3{0,OPe{0,OdO, (10) 

where a and /3 design the cartesian components. By def- 
inition, we have tr{F{^)) = 1. Introducing the harmonic 
expression (|9]) in (fTO)) . we get 

ac(0 =2(Fi(C)-F2(O)cos2[0,(e)-0.], (11) 

where the subscripts 1 and 2 refer to the principal values 
of F{S^) and 0c{£,) represents the privileged direction of 
the partial fabric tensors F{£^). Note that, up to statis- 
tical fiuctuations, the principal directions of the fabric 
and stress tensors coincide in the critical state for each 
^-network, so that the phase factor cos2[9c{£,) — 9a] is 
either equal to 1 when Od^) — 9^ or equal to —1 when 
e{0 =0. + 7r/2. 



FIG. 18: Partial fabric anisotropy (Xc 3bS 3. function of force 
cutoff ^ normalized by the mean force (/) for different values 

of TJ. 

FigHH displays as a function of ^ for all values of 
rj. For the disk packings (77 = 0), the anisotropy of weak 
contacts is negative but increases in absolute value and 
reaches its peak value at ^ ^ 1. This negative value 
indicates that in disk packings the weak contacts are ori- 
entated preferentially perpendicular to the major prin- 
cipal stress direction [ll|. As more contacts come into 
play with increasing ^, the partial anisotropy ac(^) be- 
comes less negative and finally changes sign, showing that 
the strong contacts are mainly along the major prin- 
cipal stress direction. This bimodal behavior of stress 
transmission is a nontrivial organization of the force net- 
work and holds also in 3D in the case of sphere pack- 
ings [1^. However, it is remarkable that for elongated 
particles (ry > 0), the partial anisotropics of both weak 
and strong networks are positive, as observed in Fig. 1181 
This means that, in contrast to the disk packings, the 
weak and strong contacts in packings of elongated parti- 
cles can not be differentiated on the basis of their roles 
in the ^-networks. Physically, this behavior may be in- 
terpreted by stating that the static equilibrium of the 
chains of elongated particles does not require the stabi- 
lizing effect of the weak contacts. A similar result was 
observed by Estrada et al. for disk packings at large val- 
ues of rolling resistance, which allows for the equilibrium 
of long chains of particles inter-connected by only two 
contacts [T3|. But, as we shall see below, for our elon- 
gated particles the differentiation between the two net- 
works operates via the forces carried by the ^-networks. 

The information involved in the angular distribution 
Pe may be enriched by accounting for the branch vectors 
£ which, as seen in section IIIIl refiects both the parti- 
cle size distribution and local contact modes. We thus 
consider here the average normal and tangential branch 
vector components {in){9,0 and (^t)(0,O defined in 
obtained by averaging and it over the contacts ori- 
ented along 9 within a centered angular interval A^. As 
for Pe, we evaluate these functions in the critical state, 
for different values of rj and as Fig. [12] shows the 
functions (in) {9,0 ^^'^ {^t){9,0 i^i polar coordinates for 
77 = 0.5 and for several values of ^. These functions are 
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anisotropic with an anisotropy which depends on ^. We 
introduce here their truncated expansion on an orthonor- 
mal Fourier basis: 



(a) 



(12) 



where a;„(^) and ait{S,) are the normal and tangential 
branch anisotropics in the ^-networks. Note that, by 
construction we have an = for disks (77 = 0). The 
analytical form of (it) {9, £,) results from the orthonormal 
nature of the Fourier basis and the fact that the mean 
value of it vanishes due to disorder: 



{et){0,O Pe{0,0 d9 = 0. 



(13) 



Fig. [19] shows that this functional form provides a good 
approximation of the data. 

For the calculation of ai„(£) and ait(^), we introduce 
the following branch tensors [22| : 



X%{0 = ^J{it){O,Onc.{Ot0{OPe{0,Od9, 



The following relations are then easily obtained: 



(14) 



where — x^' 



X 



-X^"(0]/tr-[x'"(^)]-«c(0, 

x^(0]/t'^[x'M]-«c(0-«;n(0, 

(15) 

and the subscripts 1 and 2 refer 



to the principal values of each tensor. By construction, 
we have trx^ = (Xi + X2) = (^)- Note also that the two 
partial branch vector anisotropies a/„ and an may take 
positive or negative values depending on the orientations 
9in and 9it of the two tensors with respect to 9a. 

Figure BUI shows the branch- vector anisotropies a/„(^) 
and ait{£,) as a function of ^ in the critical state for all 
values of 77. a/„(^) is positive for 77 = and 77 = 0.1 and 
increases slightly with ^, but for more elongated particles 
it takes negative values, which means that the particles 
tend to form longer branch vectors with their neighbors in 
the direction of extension. As f increases, this anisotropy 
increases in absolute value and reaches a plateau after 
passing by a peak value at a point in the range ^ G [1, 2]. 
This behavior suggests that the particles touch preferen- 
tially along their minor axes when the contact orientation 
is close to the compression axis (in the strong network), 
and along their major axis when the contact orientation 
is close to the extension axis (in the weak network), in 
agreement with the fact that the longest branches are in 
the weak network ; see Sec. IIIll As to a;t(^), its value 
is always negative and increases monotonically with ^ in 
absolute value. Note also that, for all values of ait{£,) 
is much higher than a/„ (^) while both remain weak com- 
pared to adO- 





FIG. 19: Distributions of (^„)(6l,^)(a) and (^t>(6',C)(b) (sym- 
bols) in polar coordinates for r; = 0.5 and several values of the 
force cutoff ^ together with their Fourier fits (|19p (full lines). 




FIG. 20: Partial normal and tangential branch vector length 
anisotropies aj„ and au as a function of force cutoff ^ normal- 
ized by the mean force (/) for different values of 77. 



B. Force anisotropies 

We now consider the angle-averaged normal and tan- 
gential forces, {fn){d,0 ^i^d {ft}{&,£.) in the ^-network. 
A second order Fourier expansion provides an adequate 
representation of these distributions for all values of ^ as 
shown in Fig. [5TJ 



{fn)iO,0 = (/){l + a/„(0 cos 2(0-0,)} 
{ft){0,O = {f)aftiOsin2{9-9,), 



(16) 



where a/ri(C) ^-nd a/t(C) ^'^^ the amplitudes of normal 
and tangential force anisotropies in the ^-networks. No- 
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FIG. 21: Distributions of {U){0,^)(a) and (/t)(6', C)(b) (sym- 
bols) in polar coordinates for 77 = 0.5 and several values of the 
force cutoff ^ together with their Fourier fits (|16|1 (full lines). 



tice that we have (/t) = as a consequence of weak cor- 
relation between the branch vectors and contact forces 
as shown in Fig. [13] and the balance of force moments. 
Morevover, the orthogonality between the normal and 
tangential forces implies that the peak value of (/t)(6',C) 
occurs at an angle rotated by 7r/4 with respect to that of 
are rotated to those of (/n)(6',C)- 

As for the branch length vectors, the calculation of the 
anisotropy parameters a/n(C) and a ft (£) can be done by 
means of the following force tensors pH. [23|: 

(17) 

With these definitions, the following relationships can 
easily be established: 



tr[x-'"(w)] 
tr[x^(cx))] 



(18) 

adO-afniO, (19) 



where = X^"' + X'^* s-nd the indices 1 and 2 refer to 
the principal values of each tensor. By construction, we 

have tr{x-^) = Xi + X2 — if)- The two partial force 
anisotropics a/„ and aft may take positive or negative 
values depending on the orientations 9fn and Oft of the 
two tensors with respect to Oa-- 

The normal and tangential force anisotropics are plot- 
ted in Figl^H as a function of ^ for all values of 77. A 
remarkable feature of a fn {£,) is that its value is negative 
in the weak network (^ < 1) for all elongated particles, 
i.e. for all values of ry in exception to 77 = where it 
remains positive for all ^. Hence, the weak forces in a 
packing of elongated particles occur at contacts prefer- 
entially oriented orthogonally to the principal stress di- 
rection 9^ whereas in a disk packing they are parallel. 
As we saw before, an inverse behavior occurs for the con- 
tact anisotropics, i.e. the weak contacts in the packings 
of elongated particles are parallel to the principal stress 
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FIG. 22; Partial normal and tangential force anisotropics a/„ 
and aft as a function of force cutoff ^ normalized by the mean 
force {/) for different values of 77. 



direction and orthogonal for the disk packings. a/„(^) 
increases in absolute value as ^ increases and passes by 
a peak at exactly ^ = 1, then declines as more contacts 
from the strong network with a positive contribution to 
the anisotropy are included in the ^-network. At larger 
values (beyond ^ ~ 2 for nearly all values of 77), a/„(^) 
becomes positive as the strong forces tend to be paral- 
lel to the principal stress direction. This unmonotonic 
behavior of the partial force anisotropics for the elon- 
gated particles and the partial contact anisotropics for 
the disk packings underlies the differentiation between 
the weak and strong networks according to the values of 
the normal contact forces with respect to the mean force 
(^ = 1). The difference between the elongated particle 
packings and disk packings reflects the formation of side- 
side contacts oriented along the principal stress direction 
tending to capture the strong force chains. 

The tangential force anisotropy a/t(^) is an increasing 
function of both ^ and 77. Its value is generally below 
a/„(^), but becomes comparable for the most elongated 
particles for which the friction mobilization plays a key 
role as discussed previously. This is plausible as the tan- 
gential force anisotropy represents friction mobilization 
at contacts oriented at 7r/4 with respect to the major 
principal stress direction. 
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C. Stress tensor 

The physical importance of geometrical and mechan- 
ical anisotropies becomes clear when it is considered in 
connection with the stress tensor. As shown by Eq. |31 the 
stress tensor is a function of discrete microscopic parame- 
ters attached to the contact network. It is also possible to 
attribute a stress tensor to each ^-network by restricting 
the summation to the corresponding contacts: 

't(0-^E/"(^)^M^)- (20) 

cev 

For sufficiently large systems, the dependence of volume 
averages on individual discrete parameters vanishes [22l . 
l45j and the discrete sums can be replaced by integrals as 
follows: 

<Jam ^nJ UiOe^iO PifiOdf d£, (21) 
Jo 

where Pif is the joint probability density of forces and 
branch vectors in the ^-networks, Uc is the number den- 
sity of contacts for the whole system and is the inte- 
gration domain in the space {£, /). 

The integral appearing in Eq. ((2T|) can be reduced by 
integrating first with respect to the forces and branch 
vector lengths. Considering the components of the forces 
and branch vectors in contact frames (n, t), and neglect- 
ing the branch/force correlations (see Fig lTB)) . we get 

TT 



X P{e,£.) de. (22) 

The expression of the stress tensor by this equation makes 
appear explicitly the average directional functions repre- 
senting the fabric and force states. 

Using the harmonic approximation introduced before, 
Eq. ((22|) can be integrated with respect to space direction 
9 and we get the following simple relation: 

^ ^ UadO + ainiO + MO + cifniO+aMO}, (23) 
P ^ 

where the cross products among the anisotropy param- 
eters have been neglected. This relation expresses the 
normalized shear stress as a half-sum of texture and 
force anisotropies. Fig l23] displays the partial shear stress 
q{Ct/p s-s a function of ^ together with the approximation 
given by Eq. [531 As we see, equation (^5)) provides an 
excellent fit to the data for all values of ^ and 77. Interest- 
ingly, g(^ < V)/p is zero for disk packings, implying that 
strong forces carry the whole deviatoric load. The partial 
stress deviator g(^ = l)/p in the weak network increases 
slightly with 77 but remains in all cases weak (below 0.1). 




FIG. 23: Partial shear stress q/p as a function of force cutoff 
5 for different values of t] (plain line) together with approxi- 
mation given by Eq. [23] (points). 

This transition reflects a qualitative change in the condi- 
tion of local force balance in the presence of clusters as 
shown in Fig. [TTJ In other words, for these packings the 
weak network sustains also partially the deviatoric load 
applied to the system. The weak values of q/p in the 
weak network is a consequence of the large positive value 
ac{S, = 1) — 0.3 which compensates the negative values 
of a/„($ = 1), a;„(^ = 1) and ait{£, = 1). 

V. SUMMARY 

In summary, using contacts dynamics simulations, we 
analyzed the granular texture and topology of forces 
chains in various packings composed of elongated par- 
ticles under biaxial compression. As compared to disk 
packings, the effect of particle elongation is to enhance 
the heterogeneity of the packings by the clustering of the 
particles according to their contact modes. In particu- 
lar, the side/side contacts tend to capture strong force 
chains and be oriented orthogonally to the major princi- 
pal stress direction. These features are reinforced as the 
particle elongation is increased. The probability densities 
of the normal forces become broader with stronger force 
chains characterized by an exponential distribution as in 
disks packings, and with higher number of weak forces 
decreasing as a power law with the force. 

An interesting finding of this work concerns the dif- 
ferentiation between the strong and weak force networks 
for elongated particles. In contrast to disks packings, 
where the contacts in the weak network are on the aver- 
age perpendicular to the contacts in the strong network, 
the contacts in a packing of elongated particles are, on 
the average, oriented along the major principal stress di- 
rection both in the weak and strong networks. But, the 
weak forces in the case of elongated particles show a neg- 
ative anisotropy in the sense that the average normal 
force in the weak network has its maximum value in the 
contacts perpendicular to the strong network. In other 
words, while in the disk packings the strong forces chains 
are propped by many weak lateral contact, for elongated 
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particles the strong force chains are laterally sustain by 
less contact but larger weak forces. A harmonic decom- 
position of the stress tensor shows, however, that for both 
disks and elongated particles, the compensating effects of 
force and contact anisotropics lead to small shear stress 
deviator carried by the weak network. 

Our simulation data indicate that the larger global 
shear strength of a packing of elongated particles in- 
creases with elongation mainly due to the increase of 
friction mobilization and friction force anisotropy. The 
normal force anisotropy is large but nearly independent 
of elongation. On the other hand, the correlation be- 
tween contact forces and branch vectors joining particle 
centers reveal a sub-network of weak contacts with hight 
friction mobilization and small branch vector length. 

In conclusion, the packings of elongated particles in 
2D reveal a nontrivial texture allying the geometry of 
the particles with the preferred orientations of the con- 
tacts induced by shearing and equilibrium of particles. 



Some features are reminiscent of disk packings but are 
strongly modulated by the particle shape. More work is 
underway to clarify the effect of particle shape by focus- 
ing on the local structures. On the other hand, many 
aspects of the packings analyzed in this paper are spe- 
cific to two dimensions. The side/side contacts in 3D be- 
tween particles of spherocylindrical shape do not give rise 
to nematic ordering and the particle rotations and forces 
moments play a major role in the equilibrium of such 
particles. This point can only be analyzed by performing 
3D simulations of large packings of sphero-cylinders of 
varying elongation. However, since the class of side/side 
contacts controls to a large extent the specific behavior 
of elongated particles in 2D, we believe that similar fea- 
tures should occur in 3D for platy particles, which may 
give spontaneously rise to geometrical chains of face/face 
contacts. Such simulations require, however, much more 
computational effort. 
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